ON THE STABILITY OF STOCHASTIC JUMP KINETICS 



STEFAN ENGBLOM 

Abstract. Motivated by the lack of a suitable framework for analyzing pop- 
ular stochastic models of Systems Biology, we devise conditions for existence 
and uniqueness of solutions to certain jump stochastic differential equations 
(jump SDEs). Working from simple examples we find reasonable and explicit 
assumptions on the driving coefficients for the SDE representation to make 
sense. By 'reasonable' we mean that stronger assumptions generally do not 
hold for systems of practical interest. In particular, we argue against the 
traditional use of global Lipschitz conditions and certain common growth re- 
strictions. By 'explicit', finally, we like to highlight the fact that the various 
constants occurring among our assumptions all can be determined once the 
model is fixed. 

We show how basic perturbation results can be derived in this setting such 
that these can readily be compared with the corresponding estimates from 
deterministic dynamics. The main complication is that the natural path-wise 
representation is generated by a counting measure with an intensity that de- 
pends nonlinearly on the state. 



1. Introduction 

The observation that detailed modehng of biochemical processes inside living 
cells is a close to hopeless task is a strong argument in favor of stochastic models. 
Such models are often thought to be more accurate than conventional rate-diffusion 
laws, yet remain more manageable than, say, descriptions formed at the level of in- 
dividual molecules. Indeed, several studies [19, 27, 28] have showed that noisy 
models have the ability to capture relevant phenomena and to explain actual, ob- 
served dynamics. 

In this work we shall consider some 'flow' properties of a stochastic dynamical 
system in the form of a quite general continuous-time Markov chain. Since the 
pioneering work of Gillespie [12, 13], in the Systems Biology context this type of 
model is traditionally described in terms of a (chemical) master equation (CME). 
This is the forward Kolmogorov equation of a certain jump stochastic differential 
equation (jump SDE for brevity), driven by independent point processes with state- 
dependent intensities. Despite the popularity of the master equation approach, little 
analysis on a per trajectory-basis of actual models has been attempted. 

In the general literature, when discussing existence/uniqueness and various types 
of perturbation results, different choices of assumptions with different trade-offs 
have been made. One finds that the treatment often falls into one of two categories 
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taking either a "mathematical" or a "physical" viewpoint. Either the conditions are 
highly general but with subsequently less transparent proofs and resulting in more 
abstract bounds. Or the conditions are formed out of convenience, say, involving 
global Lipschitz constants, and classical arguments carry through with only minor 
modifications. 

Protter [26, Chap. V] offers a nice discussion from the mathematical point of view 
and in ascending order of generality, including the arguably highly unrestrictive 
assumption of locally Lipschitz continuous coefficients. Other authors [2, Chap. 6], 
[29, Chap. 3-5] also treat the evolution of general jump-diffusion SDEs in continuous 
state spaces. 

A study of the flow properties of jump SDEs is found in [23], where the setting is 
scalar and the state continuous. In [14] jump stochastic partial differential equations 
are treated, and existence/uniqueness results as well as ergodic results for the case 
of a multiplicative noise, are found in [22, 21]. Numerical aspects in a similar setting 
are discussed in [11]. 

In a more applied context, stability is often thought of as implied from physical 
premises and the solution is tactically assumed to be confined inside some bounded 
region [17, Chap. V]. The fundamental issue here is that for open systems in a 
stochastic setting, there is a non-zero probability of reaching any finite state and 
global assumptions must be formed with great care. The analysis of open networks 
under an a priori assumption of boundedness is therefore quite difficult to inter- 
pret other than in a qualitative sense. Notable examples in this setting include 
time discretization strategies [15, 20], time-parallel simulation techniques [7], and 
parameter perturbations [1, manuscript]. 

Evidently, essentially no systems of interest satisfy global Lipschitz assumptions 
since the fundamental interaction almost always takes the form of a quadratic 
term. Interestingly, for ordinary differential equations, it has been shown [18] 
that Lipschitz continuous coefficients imply a computationally polynomial-space 
complete solution; thus providing a kind of explanation for the convenience with 
this weak feedback assumption. It is also known [16], that with SDEs, superlincarly 
growing coefficients may in fact cause the forward Euler method to diverge. 

1.1. Agenda. The purpose of this paper is to devise simple conditions that imply 
stability for finite and, in certain cases, infinite times, and that, when applied to 
systems of practical interest, yield explicit expressions for the associated stability 
estimates. As a result the framework developed herein applies in a constructive way 
to any chemical network, of arbitrary size and topology, formed by any combination 
of the elementary reactions (2.4)-(2.7) to be presented in Section 2. Additionally, it 
will be clear how to encompass also other types of nonlinear reactions that typically 
result from adiabatic simplifications. 

As an argument in favor of this bottom- up approach one can note that, for 
evolutionary reasons, biochemical systems tend to operate close to critical points 
in phase-space where the efficiency is the highest. As we shall demonstrate in an 
explicit example to be worked out in Section 4.5, for such dynamical systems, an 
analysis by analogy might be highly misleading. 

We also like to argue that our results are of interest from the modeling point 
of view. Due to the type of phenomenological arguments often involved, judging 
the relative effect of the (non-probabilistic) epistemic uncertainty is a fundamental 
issue which has so far not rendered a consistent analysis. 
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1.2. Outline. Section 2 is devoted to formulating the type of processes we are 
interested in. We state the master equation as well as the corresponding jump 
SDE and we also look at some simple, yet informative actual examples. Since it is 
expected that the flow properties of the stochastic dynamics are somehow similar to 
those of the deterministic version, we search for a set of minimal assumptions in the 
latter setting in Section 3. Techniques for finding explicit values of the constants 
occurring among our assumptions are also devised. The main results of the paper 
are found in Section 4 where we put our theory together and prove existence and 
uniqueness results, as well as some basic perturbation estimates. The findings are 
finally applied to a qualitative analysis of a nonlinear model which displays an 
unusual sensitivity to parameter perturbations. A concluding discussion is found 
in Section 5. 

2. Stochastic jump kinetics 

In this section we start with the physicist's traditional viewpoint of pure jump 
processes and write down the governing master equations. These are evolution 
equations for the probability densities of continuous-time Markov chains over a dis- 
crete state space. Although the application considered here is mesoscopic chemical 
kinetics, identical or very similar stochastic models are also used in Epidemiology 
[5], Genetics [10] and Sociodynamics [8], to name just a few. 

We then proceed with discussing a path-wise representation in terms of a sto- 
chastic jump differential equation. The reason the sample path representation is 
interesting is the possibility to reason about flow properties and thus compare func- 
tionals of single trajectories. This is generally not possible with the master equation 
approach. 

For later use we conclude the section by looking at some prototypical models. 
A simple analysis shows, somewhat surprisingly, that an innocent-looking example 
produces second moments that grow indefinitely. 

2.1. Reaction networks and the master equation. We consider a chemical 
network consisting of D different chemical species interacting according to R pre- 
scribed reaction pathways. At any given time t, the state of the system is an integer 
vector X{t) G = {0,1,2,...}-° counting the number of individual molecules 
of each species. A reaction law is a prescribed change of state with an intensity 
defined by a reaction propensity, Wr ■ Z;^ — ^ R_|_ . This is the transition probability 
per unit of time for moving from the state a; to a; — N,.; 

(2.1) P [Reaction r occurs in [0, dt]\ X{t) ~ x] — Wr{x) dt + o{dt). 

We specify such a reaction by the convenient notation 



where G Z is the transition step and is the rth column in the stoichiometric 
matrix N G Z^^^. Informally, for states x{t) G R;^, we can picture (2.1)-(2.2) as 
a stochastic version of the time-homogeneous ordinary differential equation 

R 

(2.3) x'{t) = - ^ N^u;^(a;) -nw{x) F{x), 

r=l 

where w{x) = [wi{x), . . . ,ijUfi{x)]'^ is the column vector of reaction propensities. 
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The physical premises leading to a description in the form of discrete transition 
laws (2.2) often imply the existence of a system size V (e.g. physical volume or 
total number of individuals). For instance, in a given volume V the elementary 
chemical reactions can be written using the state vector x = [a, b]'^ , 



(2.4) ^ A, 

(2.5) A 

(2.6) A + A ''^"('--^^/^ 0, 

(2.7) A + i3^^0, 



with the names of the species in capitals. These propensities are generally scaled 
such that Wrix) — Vur{x/V) for some dimensionless function Ur- Intensities of 
this form are called density dependent and arise naturally in a number of situations 
[9, Chap. 11]. For the rest of this paper, we conveniently take V =^ 1 and defer 
system's size analysis to another occasion. 

The models we consider here all have states in the positive integer lattice and 
the assumption that no transition can yield a state outside Zi^ is therefore natural. 
We make this formal as follows [4, Chap. 8.2.2, Definition 2.4]: 

Assumption 2.1 {Conservation and stability). For all propensities, Wrix) = 
for any x S such that x — Nr ^ Z^, and we also restrict initial data to Z;^. 
Furthermore, for all finite arguments x, we assume that w^^x) is finite. 

To state the chemical master equation (CME), let for brevity p{x, t) = P(X(t) = 
x\ X(0) = So) be the probability that a certain number x of molecules is present 
at time t conditioned upon an initial state X{{)) — xq. The CME is then given by 
[17, Chap. V] 

(2.8) ^ ^ + n^)p{x + Nr, t) - Wr{x)p{x, t) =: yfp{x, t). 

The convention of the transpose of the operator to the right of (2.8) is the stan- 
dard mathematical formulation of Kolmogorov's forward differential system [4, 
Chap. 8.3] in terms of which M is the infinitesimal generator of the associated 
Markov process. This is also the adjoint of the master operator in the sense 
that {Ml^p, q) ~ (p, Mq) in the Euclidean inner product over the state space. An 
explicit representation is 

R 

(2.9) M(7 = ^u;^(a;)[g(2:-N,.)-q(2:)]. 

r=l 

Under assumptions to be prescribed later on it follows that the dynamics of the 
expected value of some time-independent unknown function /, conditioned upon 
the initial state Xq, can be written 

|£;-[/(A,)]= 5: ^P^f{x)^{M^p,f)^ 

R. 

(2.10) = (p, Mf ) = ^ [WriXt) {f{Xt - Nr) - f(Xt))] . 

r=l 
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We now consider a path- wise representation for the stochastic process Xf. 

2.2. The sample path representation. In the context of analyzing models in 
stochastic chemical kinetics, the path-wise jump SDE representation seems to have 
been first put to use in [25, unpublished manuscript], and it was later further detailed 
in [20]. We thus assume the existence of a probability space P) with the 

filtration J^t>o containing i?-dimensional Poisson processes. The state of the system 
X{t) G will be constructed from a stochastic integral with respect to suitably 
chosen Poisson random measures. 

The transition probability (2.2) defines a counting process TTr{t) counting at time 
t the number of reactions of type r that has occurred since t — 0. It follows that 
these processes fully determine the state X{t), 

B. 

(2.11) Xt=Xo-Y,Nr7Tr{t). 

The counting processes are obtained from the transition intensities (cf. (2.1)) 

(2.12) P[Trr{t + dt) - TTr{t) = 1| Ft] = Wr{Xt-)dt + o{dt) , 

where by X{t—) we mean the value of the process prior to any transitions occurring 
at time t, and where the little-o notation is understood uniformly with respect to 
the state variable. Alternatively, using a random time change representation [9, 
Chap. 6.2], we can produce the counting process from a standard unit-rate Poisson 
process 11^, 

(2.13) TTr{t) = n,. Wr{Xs-) ds^ . 

The marked counting measure [3, Chap. VIII] ^r{dt x dz; cj) with w € defines 
an increasing sequence of arrival times G R+ with corresponding "marks" Zi S 
/ [0, 1] according to some probability distribution which we will take to be 
uniform. The intensity mj.{dt x dz) of fi^idt x dz) is the Lebesgue measure scaled by 
the corresponding propensity, mr{dt x dz) = Wr{Xt-) dt x dz. Using this formalism, 
(2.11) and (2.13) can be written in the jump SDE form 

(2.14) dXt = - J^nfj,{dtxdz), 

where /j, — [fii, . . . , hb]'^ ■ Here, the time t ~ t to the arrival of the next reaction 
of type r is exponentially distributed with intensity Wr{Xt-). Note that, by virtue 
of the nature of the propensities, the intensities of the counting processes therefore 
depends nonlinearly on the state [3, Chap. II. 3]. 

Using that the point processes are independent and therefore have no common 
jump times [4, Chap. 8.1.3] , we can obtain a more transparent notation in terms of a 
scalar counting measure. Define for this purpose and for any state x the cumulative 
intensities 

r 

(2.15) Wr{x) = Y,^s{x), 

s=l 

such that the total intensity is given by W{x) = Wji{x). Let the marks Zi be 
uniformly distributed on I. Then the frequency of each reaction can be controlled 
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through a set of mdicator functions Wr : x / — > {0,1} defined according to 

[Z.W) Wr(x, - I Otherwise. 

Put w{x) = [wi(a;; z), . . . , w/j(a;; z)]-^ and define also for later use the indicator 
form 

(2.17) F{x; z) = -Nw{x; z), 
such that 

(2.18) F{x) = ^ F{x; z)W{x) dz, 

where F{x) is defined in (2.3). 

The jump SDE (2.14) can now be written in terms of a scalar counting random 
measure through a state-dependent thinning procedure [6, Chap. 7.5], 



(2.19) dXt^- J NwiXt-; z)fi{dtxdz). 

Eq. (2.19) expresses exponentially distributed reaction times that arrive according 
to a point process of intensity m{dt x dz) = W{Xt~) dt x dz carrying a mark which 
is uniformly distributed in /. This mark implies the ignition of one of the reaction 
channels according to the acceptance-rejection rule (2.16). 

We will frequently decompose (2.19) into its "drift" and "jump" parts, 

(2.20) dXt = -Nw{Xt)dt- J^Nw{Xt-; z){ji~m){dt x dz). 

The second term in (2.20) is driven by the compensated measure (/i — m) and is a 
local martingale provided in essence that the path is absolutely integrable (see [3, 
Chap. VIII. 1, Corollary C4] for details). 

2.2.1. Localization; Ito's and Dynkin's formulas. In analytic work it is often nec- 
essary to 'tame' the process by deriving results under a stopping time rp := 
inft>o{|l^t|| > P} in some norm. Results for the stopped process Xtf\Tp can then 
be transferred to the original process by letting P — > oo under suitable conditions. 

Although there are many general versions of Ito's change of variables formula 
available in the setting of semi-martingales (see for example [29, Chap. 2.7] and [26, 
Chap. 11.7]), we shall get around with the following simple version [2, Chap. 4.4.2]. 
By the properties of the semi-martingale pure jump process we have for t = t Arp 

(2.21) 

f{Xi) - f{Xo) = /(^^) - /(^-) = / / /(^.^) - f(^s-) Kds X dz), 

Q<s<i " ^ 

where the sum is over jump times s G (0,i]. Using that Xs = Xs- — Nw{Xs-. z) 
we can write this in differential form as 

(2.22) dfiXt) = f f{Xt_ - nw{Xt_- z)) - f{Xt_)^i{dt X dz). 
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Alternatively, decomposing (2.21) into drift- and jump parts and taking expectation 
values we get, since the compensated measure is a local martingale, 

Ef{Xi) - EfiXo) f{Xs- - m{Xs-; z)) - f{X,^)m{ds x dz) 

= E j J [fiX,_ -NwiXs-; z)) - fiXs-)]WiX,_)ds X dz 

(2.23) = ^ / E [(-^(^^ - - f{Xs))wr{X,)] ds. 

This is Dynkin's formula [4, Chap. 9.2.2] for the stopped process and we note that 
(2.10) is just a differential version. 

2.2.2. The validity of the master equation. With this much formalism developed, 
we may conveniently quote the following result: 

Theorem 2.1 {[4, Chap. 8.3.2, Theorem 3.3]). Under Assumption 2.1, and if 
additionally, for t € [0, T] it holds that 

(2.24) EW{Xt) <oo, 

then (2.8) is valid for t G [0,r]. 

Since the governing equation (2.10) for the expected value of f{Xt) is a direct 
consequence of (2.8), we can similarly conclude the following: 

Corollary 2.2. Under the assumptions of Theorem 2.1, and if, moreover, in an 
arbitrary norm \\ • |j, 

(2.25) E\\f{Xt)\\<^, te[o,r], 

then (2.10) is valid for t e [0,T]. 

In stating these results we have suppressed the conditional dependency on the 
initial state which we for simplicity consider to be some non-random state xq. 

2.3. Concrete examples. Consider the bi-molecular birth-death system, 

(2.26) S L 

that is, the system is in contact with a large reservoir such that A- and B-molecules 
are emitted at a constant rate fci. Additionally, a decay reaction happens with 
probability fc2 per unit of time whenever two molecules meet. For this example we 
have the stoichiometric matrix 



and the vector propensity function 

w{x) — 

where x — [xi,X2]'^ — [a,b]'^. 



[fci,fci,fc2a6]^. 
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For a state Xt = [At.Bt]^, define Ut ^ At - Bt ^ Z. Ito's formula (2.22) with 
f{x) ~ xi — X2 yields 



(2.27) 



dUt = dfiXt) 



[-1, l,0]w{Xt-; z) fi{dt X dz), 
which upon a moments consideration is just the same thing as the model 
(2.28) 



fcl 



that is, a constant intensity discrete random walk process. An explicit solution is 
the difference between two independent Poisson distributions, 

(2.29) Ut = UQ + Ui{kit)-U2{kit) r^M{Un,2kit), ast^oo, 

where A/" is a normally distributed random variable of the indicated mean and 
variance. Hence Ut fluctuates between arbitrarily large and small values as t ^ oo. 

2.3.1. Reversible versions. From time to time below we shall be concerned with the 
following closed version of (2.26), consisting of a single reversible reaction. 



(2.30a) 



A 



kiab 

B U C 

k2C 



This is clearly a finite system since the number a 
open version of the same system is 



6 + 2c is always preserved. An 



(2.30b) 



A 



B U C ^\ 

k^C fc4 



and will prove to be a useful example in the stochastic setting since formally, all 
states in 7i\ are reachable. For (2.30a) we have 



(2.31) 

while (2.30b) is represented by 
(2.32) 



w{x) = [kiab, fee] 



1 
1 
1 1 



^(a:;) = [kiab, fee, fee, fe] 



These examples, while very simple to deal with, will provide good counterexamples 
in both Section 3 and 4. 



3. Deterministic stability 

In this section we shall be concerned with the deterministic drift part of the 
dynamics (2.20). We are interested in techniques for judging the stability of the 
time-homogeneous ODE (2.3), the so-called reaction rate equations implied by the 
rates (2.2). Stability and continuity with respect to initial data are considered 
in Sections 3.1 and 3.2, and techniques for explicitly obtaining all our postulated 
constants are discussed in Section 3.3. 

Initially we will consider states x G R^, but we will soon find it productive to 
restrict the treatment to a; G R:^. In order to remain valid also in the discrete 
stochastic setting, however, constructed counterexamples will remain relevant also 
when restricted to Z£. 
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3.1. Stability. Many stability proofs can be thought of as comparisons with rele- 
vant linear cases. This is the motivation for the well-known Gronwall's inequality 
which we state in the following two versions. 

Lemma 3.1. Suppose that u'{t) < A + au{t) for t>Q. Then 

(3.1) u(i) < u(0)e"' + - (e"* - 1) . 

a 

The same conclusion holds irrespective of the differentiability of u but with a > 
and under the weaker integral condition 

(3.2) u{t)<u{0)+[ A + au{s)ds. 

Jo 

Proof. The constant A can be disregarded by using the substitution v{t) = u{t) + 
A/a. For the differential version one integrates the given inequality using exp(— at) 
as an integrating factor. For the integral version, define 

{v{t) - w(0)e"* < ) w{t) := v{0) + [ av{s) ds - w(0)e"*. 

Jo 

Under the assumed integral inequality we readily get w'lt) < aw{t) such that 
wit) < for t > 0. □ 

The most immediate way of comparing the growth of solutions to the ODE (2.3) 
to those of a linear ODE is to require that the norm of the driving function is 
bounded in terms of its argument; 

(3.3) \\Fix)\\<A + a\\x\\, 
since then by the triangle inequality, 

(3.4) ||a;(i)|| < llxoll + / \\F{x)\\dt<\\xo\\+ f A + a\\x{t)\\dt, 

Jo Jo 

where Gronwall's inequality applies. Unfortunately, (3.3) is a too strict requirement 
for our applications. 

Proposition 3.2. The bi-molecular birth-death system (2.26) does not satisfy (3.3). 

Proof. We compute ||i^(a;)|| ~ \\ — l^w{x)\\ = \/2\ki — k2ab\ for a state x = [a,bY . 
Hence for a = b = N — 0,1,... we have for iV large enough that ||F(a;)|| = 
^/2k2 ■ — V^ki, which can clearly never be bounded linearly in = y/2N . □ 

The problem with the simple condition (3.3) is that it does not take the direc- 
tion of growth into account; the offending quadratic propensity in (2.26) actually 
decreases the number of molecules. To deal with this, let x G be an arbitrary 
vector defining an "outward" direction. The length of the composant of the driving 
function along this direction is (x, F{x)) and in order not to have x driven too 
strongly out along this ray we may, in view of Gronwall's inequality, naturally re- 
quire that (a;/||x||, F(a:)) < constant x for ||x|| sufficiently large. Equivalently, 
for any x, 

(3.5) {x,F{x))<A + a\\x\\^, 
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from which one deduces 

(3.6) ^^^^ = {x,Fix))<A + a\\xr, 

where Gronwall's inequahty applies anew. The assumption (3.5) is weaker than 
(3.3) since the former imphes the latter by the Cauchy-Schwarz inequality. Indeed, 
as in the proof of Proposition 3.2 it is readily checked that for the bi-molecular 
birth-death system (2.26), we get {x,F{x)) = ki{a + b) — k2{a + b)ab which this 
time readily can be bounded linearly in terms of = + 5^. 

Unfortunately, in the case of an infinite state space and strong dependencies 
between the species the assumption (3.5) is also often unrealistic. 

Proposition 3.3. Neither (2.30a) nor (2.30b) admits a bound of the kind (3.5). 

Proof. As in the proof of Proposition 3.2 we look at a ray x'^ = (a, b, c) = {N, N, 3N) 
parametrized by a non-negative integer A^. For (2.30a) we compute {x,F(x)) = 
(x, —Nw{x)) — (a + 6 — c)(fc2C — kiab) = kiN^ — 3^2^^, which clearly cannot be 
bounded linearly in — llA^. The same argument applies also to (2.30b). □ 

This negative result can perhaps best be appreciated as a kind of loss of infor- 
mation about the dependencies between the species in the functional form of the 
condition (3.5). The number of A- and i?-molecules is strongly correlated with 
the number of C-molecules such that, in fact, in (2.30a) a -I- 6 -|- 2c is a preserved 
quantity. By contrast, in (3.5) the growth of ||a;|p is estimated from the sum of the 
growth of the individual elements of x as if they where independent. 

A way around this limitation can be found provided that we leave the general 
case x e R^. We therefore specify the discussion to the positive quadrant x G R;^ 
and assume from now on that it can be shown a priori that the initial data xq 
belongs to this set and that the subsequent trajectory x{t) never leaves it (compare 
Assumption 2.1). 

It then follows that ||a;t||i = {l,xt) = I'^xt, where 1 is the vector of length D 
containing all ones. This vector also defines a suitable "outward" vector for states 
X € R:^ since solutions to the ODE (2.3) cannot grow without simultaneously 
growing also in the direction of 1. 

Again, in view of Gronwall's inequality Lemma 3.1, we tentatively require that 
{l,F{x)) < constant x ||x||i for sufficiently large. Equivalently, for any x, 

(3.7) {l,Fix))<A + a\\x\\i, 
implying the bounded dynamics 

(3.8) j^\\x\\,^il,F{x))<A + a\\x\\i. 

We remark in passing that the criterion (3.7) is sharp in the sense that if the 
reversed inequality can be shown to be true, then the growth of solutions can be 
estimated from below. 

Example 3.1. As a point in favor of this approach we compute for the bi-molecular 
birth-death system (2.26), (1, F{x)) = (1, -Nw(x)) = 2fci - 2k2ab which evidently 
falls under the assumption (3.7) with {A, a) = (2fci,0). For the reversible case 
(2.30a) we similarly get {l,F{x)) — —kiab + k2C such that (3.7) applies with 
{A, a) = (0,^2). Finally, and in the same fashion, the open case (2.30b) is seen to 
be covered by letting {A, a) — {k^, (/c2 — ^3) V 0). 
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The chosen "outward" vector 1 is by no means special. Clearly, any strictly 
positive vector / may be used in its place since ||x|| i and (;) := l^x are equivalent 
norms over R;^. This is a general and useful observation as it may be used to 
discard parts of a system that are closed without any restrictions on the associated 
propensities. 

Example 3.2. For the reversible system (2.30a), we have already noted that a + 
& + 2c is a conserved quantity such that the choice / = [1, 1, 2]-^ yields d/dt\\x\\(i) = 
0. The open case (2.30b) also benefits from this weighted norm in that we get 
d/dt\\x\\(i) < 2ki. 

Example 3.3. A slightly more involved model reads as follows: 

^ A A + B 3C 1 

0^5 C "-^ % ] 

This example has been constructed such that the quadratic reaction increases ||a;||i 
and hence (3.7) does not apply. However, taking I — [3, 3, 2]^ we get 

^lkt||(o = 6fci - 2fc3C < 6fci. 

This example hints at a general technique for obtaining suitable candidates for 
the weight vector I. Simply form the matrix N2 consisting of the columns of N 
that are affected by superlinear propensities. If a vector I > annihilating these 
propensities exists, it can be found in the null-space of — N^, readily available by 
linear algebra techniques. We omit the details. 

3.2. Continuity. For well-posedness of the ODE (2.3) we also need continuity with 
respect to the initial data. We cannot ask for uniform Lipschitz continuity since 
||F(a;) — F{y)\\ < L\\x — y\\ clearly implies (3.3) which we have already refuted. For 
the same reason, a uniform one-sided Lipschitz condition {x — y,F{x) — F(y)) < 
X\\x ~ yiP cannot be assumed to hold since it implies (3.5). The problem here is 
the global nature of the estimate and it therefore seems to be reasonable to localize 
this assumption. For instance, one might ask for 

(3.9) ix-y,F{x)-F{y))<\R\\x-y\\^ whenever \\x\\, \\y\\ < R, 

presumably with some growth restrictions on A^j. Although theoretically precise, 
such an analysis is likely to be somewhat inconvenient when it comes to estimating 
actual constants in perturbation results. We shall therefore consider the following 
version, 

(3.10) ix-y,F{x)-F{y)) < {M + fi\\x + y\\,)\\x - y\\^ , 

where the form of has been restricted to better suit the present purposes. Triv- 
ially, the norms || • ||i and || • |1 are equivalent and hence the specific choice made in 
(3.10) is just a matter of convenience. Since the idea here is to use a priori bounds 
on X and y when deriving perturbation bounds, using || -Hi (or || • ||(/)) is natural. 

Theorem 3.4. Suppose that the ODE (2.3) satisfies (3.7) and (3.10) and that 
initial data Xq G R+ implies a trajectory x{t) S R;^. Then for any t S [0,T] there 
is a unique solution x(t). Moreover, define C{t; Xo^yo) = M + ^i\\xt + yt\\i, where 
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xt and yt are two trajectories associated with initial data xq, yo & R-+; respectively. 
Then 

(3.11) \\x{t) - y{t)\\ < ||a;o - yo\\ exp (^J^ C{s; Xq, yo) ds 

= \\xo - yoW [l + (M + ^(||xo + yo\\i))t + O (t^) 
Proof. Combining (3.7) with Gronwall's inequality we get the a priori estimate 

||a;(t) + yit)\\i < \\xo + yo\\i + {a\\xo + yo\\i + 2A) (e"* - l) /a. 
Hence the (bounded) solution to 

- = 2ix - y, F{x) - F{y)) < 2C{t- xo, yo)\\x - yf. 

is readily found through its integrating factor. The order estimate is a consequence 

of the fact that 

\\xs\\i -\\xo\\i ds^O{t^), 
since the trajectory is continuous. □ 







3.3. Bounds for elementary reactions. As briefly discussed by the end of Sec- 
tion 3.1, finding bounds on A and a in (3.7) as well as a suitable weight-vector 
I amounts to basic inequalities and some fairly straightforward linear algebra ma- 
nipulations. In this section we therefore consider precise bounds in (3.10) for the 
elementary propensities (2.4)-(2.7). Since (3.10) is linear in F, a reasonable ap- 
proach is to consider linear and quadratic propensities separate (constant propen- 
sities trivially satisfy (3.10) with M — ji = 0). 

Proposition 3.5 [linear case). Write a set of R linear propensities as Wr{x) = 
qjx, r = 1,...,R, each with the corresponding stoichiometric vector Nr. Then 
F{x) := -J^r^rWrix) Satisfies {x - y,F{x) - F{y)) < M\\x - yjp with M = 
A*2 [^NQ"^] in terms of the Euclidean logarithmic norm ^2 ['] and the matrix Q 
containing the vectors as columns. In particular, in the case of a single linear 
propensity and, if as is usually the case, Qj — k6jn is all-zero except for a single 
rate constant k in the nth position, then this reduces to M — k (— Nr,n + ||I^r||)/2. 

Proof. The first assertion is immediate since the smallest such constant M by def- 
inition is the logarithmic norm (see e.g. [30]). To compute /i2 [— N^g^] when q has 
the form indicated, we determine the extremal eigenvalue of — (N^gJ + <7rN^)/2. 
By the (signed) scaling invariance of the logarithmic norm we may without loss of 
generality take fc = 1. The spectral relation for an eigenpair (A, z) can be written 
as 

"^N,.jZ„ = Xzj, j = 1,2, ...,£», j 7^ n. 



2 

1 1 7. 



2 " 2 

For non-zero A the first relation can be solved for Zj . When inserted into the second 
relation, using that z„ 7^ (or otherwise z = 0), we get a quadratic equation for A 
with a single extremal root. □ 
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Example 3.4. The simple special case in Proposition 3.5 is generally sharp except 
for when there are linear reactions affecting all species considered in the model. For 
example, in a one-dimensional state space, the single decay ^ — > with propensity 
wi(a) = ka allows the optimal value M = —k. In general Z?-dimensional space, a 
chain with unit rate constants of the form Ai ^ A2 ^ ■ ■ ■ ^ Aq ^ 0, or a closed 
loop in which the last transition is replaced with Ad Ai, both admit bounds 
M < as an inspection of the Gershgorin-discs of — (N + N-^)/2 shows. 

Other than for those special examples, for the most important linear cases, Ta- 
ble 3.1 summarizes the bounds as obtained from the special case in Proposition 3.5 
(with all reaction constants normalized to unity). 

Reaction Bound on Af 


A^B (\/2-l)/2 
A^B + C (%/3-l)/2 
Table 3.1. Linear propensities and bounds of M in (3.10). 



Proposition 3.6 [quadratic case). Write a general quadratic propensity as Wr{x) = 
Sx with S a symmetric matrix. Then Fr{x) = —f^rWr{x) satisfies (3.10) with 
M = and /i = ||a;-hy||rV2 [-Nr(a; -h 2/)'^S'] < ||Nr|| US'!! . For the special case 
that Sij ~ k{dimSjn + ^jm^in)/'^ there holds < k maxjg{„j_„} (— N,r.j + ||N,.||)/4. 

Proof. Since S is symmetric we have x'^ Sx — Sy = {x + y)'^ S{x — y). Hence an 
explicit expression for fi is obtained as follows: 

fi< sup \\x + y\\^^\\x~y\\^'^{x-y,-Nr(x + y)'^S{x-y)) 

= sup sup ||w||j'^||z;||"^(w, -NrW^S't;) < sup ^2 [-NrU^S"] . 

||«||i<i 

The indicated upper bound is derived from the fact that |/i2 [•] | < |1 • || [30]. For the 
useful special case, define first the vector q = qi + q2 in terms of qij = kSjmixn + 
yn)/'2, and (72 j — kSjn{xm + ym)/2. Using the fact that the logarithmic norm is 
sub-additive we can reuse the calculation in the proof of Proposition 3.5, 

/i2 [~Nr{x + yfS] - ^12 [-Nrq'^] < M2 [-^rql] + [-^rql] 

= {xn + yn)k{~^r,m + ||N^||)/4 + (x^ + y,„)fc(-N,.,„ + ||N^||)/4. 

□ 

Example 3.5. The most important quadratic cases are summarized in Table 3.2. 
For the dimerizations in the lower half of the table there is also a linear part M in 
(3.10). 

Example 3.6. The bi-molecular birth-death model (2.26) admits the constants 
(M,^) = (0,fc2(V2- l)/4) in (3.10). Similarly, the reversible cases (2.30a) and 
(2.30b) both obeys (3.10) with (M,/i) = (%/3- l)/2 x (A:2,fci/2). AU these resuhs 
are sharp except for the open case (2.30b) for which one can obtain a slightly smaller 
constant M by using the general formula stated in Proposition 3.5. 
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Reaction 



Bound on M 



A + B ^ A + C 



A + B-^C 
A + B^ A 



(y3-l)/4 

1/4 

^/2/4 



A + A^% 
A + A^B 



2 
75/2-1 V 



75/2 + 1 



Table 3.2. Quadratic propensities and bounds of M, /x in (3.10). 



Example 3.7. As a highly prototypical example we consider the following natural 
extension of the bi-molecular birth-death model (2.26), 



where in this example it is informative to consider the dependence on the system's 
size V . It is straightforward to show the bounds [A, a) < {2kiV, — fcs) in (3.7) and 
hence that the system is effectively bounded despite being of open character. This 
is seen from the fact that, for states > 2ki/k3-V , the dynamics is dissipative in 
the II • Ill-norm. Furthermore, from Proposition 3.5 and 3.6 we get the sharp bounds 
{M,fi) < {-kiV,k2/V ■ (\/2- l)/4) in (3.10). It follows that for states {x,y} such 
that ||a; + j/||i < 9.7fci/fc2 • V^, the dynamics is contractive in the Euclidean norm. 
For density dependent propensities we expect that ||a::|| ^ V in any norm as V 
grows, and hence the region of contractivity grows in a relative sense. Intuitively 
one expects that these results offer an insight into the evolution of the process that 
is relevant also in the stochastic setting. 



We now consider the properties of the stochastic jump SDE (2.19). For con- 
venience we start by collecting all assumptions in Section 4.1. In the stochastic 
setting the requirements for existence and uniqueness are slightly stronger than in 
the deterministic case such that the one-sided bound (3.10) needs to be augmented 
with an unsigned version, implying the assumption of at most quadratically growing 
propensities. We demonstrate that this assumption is reasonable by constructing 
a model involving cubic propensities and with unbounded second moments. On 
the positive side we show in Section 4.2 that the assumptions are strong enough to 
guarantee finite moments of any order during finite time intervals. 

We prove existence and uniqueness of solutions to the jump SDE (2.19) in Sec- 
tion 4.3. A sufhcient condition for the existence of asymptotic bounds of the pth 
order moment is given in Section 4.4 where we also derive some basic perturbation 
estimates. With these in mind and as an illustration, we treat in Section 4.5 a 
model which displays a particularly dramatic response when subject to parameter 
perturbations. 

4.1. Working assumptions. We state formally the set of assumptions on the 
jump SDE (2.19) as follows: 

Assumption 4.1. For arguments x, y E Z£, 




4. Stochastic stability 
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(i) {l,F{x))<A + a\\x\\,, 

(ii) \\w{x) - w{y)\\i < (i + A||x + y\\i)\\x - y\\, 

(iii) \\w{x)\U<T + ^\\x\\l 

(iv) {x - y, F{x) - F{y)) < {M + ^l\\x + yh)\\x ~ y\\\ 

We assume the parameters {A,L,A,r,7} to be positive but allow also negative 
values of {a, Af, /x}. 

The only unique assumptions here are (i) and (ii) since the latter clearly implies 
(iii)-(iv). However, as we saw in Section 3.3, in (iv) it is often possible to find 
sharper constants M and by considering this bound in isolation. We note also that 
although (ii) is stronger than (iv) , it is in particular valid for quadratic propensities 
as can be seen from the representation used in the proof of Proposition 3.6, 

(4.1) \yj^[:,) - w^{y)\ = \{:, + y)Ts[x-y)\ < \\S\\\\x + yUx-y\\. 

4.1.1. The danger with cubic propensities. Assumption 4.1 (iii) specifies the discus- 
sion to propensities with at most quadratic growth. To show that this is natural we 
now demonstrate that care should be taken when considering cubic propensities. 

Example 4.2. Consider the model 

3X ^^^^^^^^^^^ X \ 

such that the stoichiometric vector is given by N = [—2, 1], and hence that the drift 
-Nw(x) = 0. 

Proposition 4.1. For the model in Example ^.2, if Xq > 3, then the second 
moment explodes in finite time. 

Proof. Assume that both the second and the third moment are bounded for t € 
[0, T) with T > 0. From (2.10) we get the governing equation 

j^EX^ ^E[3XtiXt-l)iXt-2)], 

such that the growth of the second moment remains bounded only provided that 
the third moment remains finite. It is convenient to look at the cumulative third 
order moment. From (2.10), 

^ECsiXt) ^EXtiXt - l){Xt ~ 2) = E [9C3{Xt){Xt ~ 2/3)]. 
at at 

By the arithmetic-geometric mean inequality, a; — 2/3 > a;— 1 > [x{x— l)(a; — 2)]^/'^, 
such that by Jensen's inequality, 

d 

We put = E C3{Xt) and get the differential inequality 

rf(lM<_3 
dt ~ 

which can be integrated and rearranged to produce the bound 

(4.2) EX,{X, - 1)(X, - 2) > 2) 



-EC3{Xt)>9E C3{Xt)*/^ >9[EC3{Xt)] 



14/3 



l-3<Xo(Xo-l)(Xo-2) 
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Hence the third, and consequently also the second moment explode for some finite 
t whenever > 3. □ 

Interestingly, we note that if Xq — 3, then the probability that the cubic decay 
transition occurs first is 1/3, and if this happens the state of the system will be 
stuck with a single molecule indefinitely. 

4.2. Moment bounds. In this section we consider general moment bounds derived 
from (2.23) using localization. To get some guidance, let us first assume that the 
differential form of Dynkin's formula (2.10) is valid. Since any trajectory {Xt)t>o 
by the basic Assumption 2.1 will belong to Z^, we may use that H-'^^tHi = {l,Xt). 
Hence from (2.10) with f{x) — {l,x) we get that 

(4.3) = (l'^(^t)) < ^ + 

by Assumption 4.1 (i). Clearly, the differential form of Gronwall's inequality in 
Lemma 3.1 applies here. A correct version of this argument unfortunately looses 
the sign of a. 

Proposition 4.2. If Assumption 4-1 (i) is true, then 

E\\Xt\\i < \\Xo\\i cMa+t) + A{exp{a+t) - l)/a+, 
where — a\/ 0. 

Here and below we shall make use of the stopping time Tp = inft>o{||Xt||i > P} 
and define t ^ t A Tp. 

Proof. From (2.23) with f{x) — (l,x) we get that 

(4.4) E\\Xi\\,^\\Xo\\i+E f il,FiXs))ds<\\Xo\\i+E f A + a+\\Xsh ds, 

Jq Jq 

By the integral form of Gronwall's inequality in Lemma 3.1 we deduce in terms of 
Yt :— Xt/\rp that 

(4.5) E\\Yt\\i < \\Xo\\i exp{a+t) + A{exp{a+t) - l)/a+ 

such that the same bound holds for Xt by letting P oo. □ 

We attempt a similar treatment for obtaining bounds in mean square. Assuming 
tactically that (2.10) is valid, writing ||a::||2 = x'^x we get after some work that 

(4.6) 1^"^*"2 = E [l^m^w{Xt) - 2X^Nw{Xt)] 

where Nf^ = (N^j)^. We expect from Gronwall's inequality that i?||Xt||^ grows at 
most exponentially with at whenever 

(4.7) l^N2w(a;) - 2x^Nw{x) <A + a\\xf. 

However, this tentative condition is often violated in practice since the second term 
—x'^Nw{x) = {x,F{x)), and since we already know from Proposition 3.3 that this 
quantity does not admit bounds in terms of even for very simple problems. 

Surprisingly, more realistic conditions arise when looking for bounds with respect 
to \\Xt\\l instead. 
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Proposition 4.3. If for some constants j3 and B, 

(4.8) {l^nfw{x) - 2l'^Nw(x)||a;||i <B + I3\\x\\l, 

(with [l^K)"^ understood elementwise) , thenE\\Xt\\\ < \\Xo\\lexp{/3^t)+B{cxp{/3+t)- 
l)//3+. 

The proof of Proposition 4.3 follows the same pattern as for Proposition 4.2, 
but using this time f{x) = \\x\\l = {l,xy in (2.23). The condition (4.8) is 
typically more realistic than (4.7) since we recognize the term — l-^Nw(a;)||a::||i = 
(1, i^(x))||a;||i, which under the evidently reasonable Assumption 4.1 (i) is < (A + 
a||x||i)||a;||i. It follows that if (l-^N)^u>(a;) grows at most quadratically with ||a;||i, 
then this assumption is sufhcient to yield bounds in mean square. Stated formally, 

Proposition 4.4. Under Assumption 4-1 (i) and (iii) the condition (4.8) of Propo- 
sition 4.3 is true with /3 = \\{l'^N)'^\\^-f + e + 2a and B = \\{l'^N)'^\\^T + e^'^A^ 
for any e > 0. 

Proof. This is straightforward: we get by the assumptions and Holder's inequality, 

(l^N)2^(x)-21^Nu;(x)||x||i < ||(l^N)2|U||z«(x)|li+2(A + a||x||i)||x||i, 
where an application of Young's inequality yields the indicated bounds. □ 

As a strong point in favor of our running assumptions we now demonstrate 
that the above reasoning can be generalized: the condition sufficient for 1st order 
stability paired with quadratic propensities implies finite time stability in any order 
moment. 

Theorem 4.5 (Moment estimate). Under Assumption 4-1 (i) o,nd (iii), and for 
any integer p > 2, 

(4.9) EWXtf, < \\Xor,exp{/3+t)+B{exp{l3+t) - 1)//?+, 

with /3+ = /3 V 0, and bounds (3 < (p - 1 + pa) + (T + 7)(2p - 2 - p)||l^N||P, and 

B < ^p + r||i^N||p„. 

Proof Using (2.23) with f{Xt) = [l'^Xt]P we get 

/■* ^ 

E\\Xi\\P^\\Xo\\1 + E TwriXs) [l^{X,~Nr)Y -[l^X^Y 
Jo ^-1 '- 



=-E I G{Xs)ds. 
Jo 



Expanding and using Young's inequality with exponent p/{p — 1) and conjugate 
exponent p, respectively, we find that 

G{X) ^yWrix)y(P)il'^x)P-^~l^Nr)^ 



r=l 3 = 1 ^J 



J=2 



<AP+{p- i+pa)\\xr, + E ( ^ ) Mr'u'^muM^ni 



<AP + ip^l+ pa)\\x\\l + (l + (2^^ - 2 - p)\\x\\1-') |ll^N|lL(r + j\\x\\i) 
<B + (3\\x\\l 
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We thus get 

EWXsK < WXn 



\P+ f B + l3+E\\Xs\\Pds, 
Jo 



where Gronwall's integral iiiequahty apphes. Letting P — > oo we obtain the result. 

□ 

As demonstrated in Examples 3.2 and 3.3, the vector 1 is not always the best 
choice for a system's "outward" direction. In such cases we note that the entire 
statement of Theorem 4.5, including the proof, holds verbatim if a suitably nor- 
malized vector Z > is used in place of 1 and the norm || • ||i is replaced with the 
equivalent norm || • ||(;). 

4.3. Existence and uniqueness. We shall now prove that the jump SDE (2.19) 
under Assumption 4.1 has a uniquely defined and locally bounded solution. To this 
end and following [29, Sect. 3.1.2], we introduce the following space of path- wise 
locally bounded processes: 

(A ^^\\ <?2 loc/yDN / ^{ti^) '■ ^ is J^t-adapted and Z^-valued such that 1 
(4.10) ^-^+^ = 1 ^sup,g[o,T] 11^(^,^)11? <ooforVT<oo j' 

Theorem 4.6 {Existence). Let Xt be a solution to (2.19) under Assumption 4.I (i) 
and (iii). Then if WX^"" < 00, {Xt}t>o e 5^^°^(Zf). 

Proof. As before we use the stopping time Tp — infi>o{||Arf ||i > P} and put 
t = t /\Tp. We get from Ito's formula 

\\Xt\\l^\\X4l + f 2{l,F{XsmXs\\i + (-l^Nf w(Xs)ds+ 

2(1,P(X,_; z))||X,_||i + il,F{Xs-; z))^ {ti-m){ds x dz). 

Since the propensities are bounded for bounded arguments (Assumption 2.1), using 
the stopping time we find that the jump part is absolutely integrable and hence a 
local martingale Mf. We estimate its quadratic variation, 

E [M]y^ = E 



E 



< E 



2(1,P(X,_; z))\\X,_\\i + (1,P(X,_; z))^'\\{ds x dz) 

[2{-l^N)\\Xs-\\i + {-1^N)^Y w{Xs-; z)WiXs-)dzds 
f |2(-l^N)||X,||i + {-l^n)^\w{X,)^'^ ds 



1/2 



<E f Ci+C2\\Xs\\lds 
Jo 



for constants Ci , C2 > depending on the network topology N and on the constants 
r and 7 in Assumption 4.1 (iii). For the drift part we have already constructed a 
suitable bound in Proposition 4.4 such that 



!^dl?<ll^oll?+ f B + p\\Xs\\lds- 
Jo 
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Taking suprenium and expectation values we get from Burkholder's inequality [26, 
Chap. IV. 4] with some constants C3, C4 > that 

E sup \\Xs\\1<\\Xq\\1+ f + sup \\X,,\\lds. 

se[0,i] Jo s'£[0,s] 

Writing = sup^^jg ^] H-^^sHi we conclude that 

E\\Xf,^,^<\\Xo\\l+ f Cs + C,E\\X\\l^,^ds. 
Jo 

By Gronwall's integral inequality we have thus shown that can be 

bounded in terms of the initial data and time t. The result now follows by let- 
ting P — > 00 and using Fatou's lemma. □ 

We note in passing that Theorem 4.6 can readily be generalized to solutions in 
the space S^^°^{Z^) for any p> I. 

Theorem 4.7 (Uniqueness). Let Assumption 4-1 hold true. Then the 

solution to (2.19) is path-wise unique. 

We shall be using the observation that, for x e Z^, we have that ||a;||i < ||a;||2 
(referred below to as the "integer inequality"). 

Proof. Under the same stopping time as before we get from Ito's formula after 
discarding the martingale part 

E\\Xi-Yif=E f 2{X,-Y,,F{Xs)-FiYs)) 
Jo 

(4.11) +{l^N^)\w{Xs)-w{Y,)\ds 

<E f 2{M + fi\\Xs+Ys\\i)\\Xs-Ysf 
Jo 

(4.12) + ||(1^N2)||^(L + \\\Xs + Y,\\i)\\Xs - Y,\\ ds. 

From the integer inequality we find that there is a constant C > depending on P 
such that 

E\\{X ~Y){tATp)f < CE\\X,-Ysfds< CE\\{X ~ Y){s A Tp)\\^ ds. 

Jo Jo 

Using that Xq = Yq and Gronwall's integral inequality we conclude that the only so- 
lution is the zero solution. Letting P — > 00 and using Fatou's lemma the statement 
is therefore proved. □ 

4.4. Basic perturbation estimates. We next aim at deriving some stability es- 
timates with respect to perturbations in initial data and of the coefhcients. To this 
effect we shall need a slightly more general Gronwall-type estimate than before. 

Lemma 4.8. Let u and a he non-negative functions in some positive interval in- 
cluding zero. Suppose that d/dtuit)"^ < a(t)u -\- /3u^ for t > 0. Then 

(4.13) u{t) < e'5/2t J^^(o) + ^J^ a(s)e-^/2,s _ 
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The same conclusion holds with /3 > and under the weaker integral condition 
(4.14) u{tf<u{Qf+[ a{s)u{s) + f3u{sf ds. 



Proof. Performing the differentiation and dividing through by u we get the imphed 
inequahty u' < a{t)/2+(3/2u. This inequaUty has an integrating factor exp(— /3/2t) 
and is therefore readily solved. To prove the second part, define 

w{t) := ( u(0)^ + / a{s)u{s) + (3u{sf ds ' 



We get under the assumptions that w'{t) < (3/2w{t) such that w{t) < 0. □ 

Note the special case of the lemma a{t) = 1, /? = 0, and u(0) = 0, implying 
the (sharp) bound u < t/2. This observation shows that the integer inequality as 
used in the proof of Theorem 4.7 is crucial for path- wise uniqueness; without it the 
integral inequality (4.12) admits growing solutions. 

Although Theorem 4.5 shows that any moments are bounded in finite time, a 
relevant question from the modeling point of view is whether the first few moments 
remain bounded indefinitely. We give a result to this effect which relies on the 
existence of solutions in S'^^°'^{Zi^) and the assumption of at most quadratic growth. 
Under these conditions the differential form (2.10) may be used (cf. Corollary 2.2) 
such that the differential Gronwall inequality applies. 

Theorem 4.9. Under Assumption 4-1 (i) CLud (iii), suppose that for an integer 
P>1, 

(4.15) 2a + 7||(l^N)2||^(p- 1) < -Kp < 0. 

Then E\\Xt\y^ is asymptotically bounded. 

Proof. We omit the case p — 1 since it follows from (4.3) under the present as- 
sumptions. The idea of the proof is to asymptotically bound E{C + ||Xt||i)^' with 
a certain positive constant C = C{p) to be decided upon below. By (2.10) we get 
withZt = C+\\Xt\\i, 

dEZf 



r=l 

= Ej2MXt)\^i^Nrpzr' + ^ p{p-mt - dA^Nry-^ 



dt 

(-l^N,) 



=1 



where Taylor's formula was used and where 9r G [0, 1] are constants. Using the 
assumptions we get the bound 



dEZf p 
* < -E 



2iA + a\\Xth)Zr' 



dt - 2 

(4.16) + ||(l^N)2|U(p - l)(r + j\\Xt\\i)iZt + ||i^N|U)f-2 
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For the first term in (4.16) we get from tlie scaled Young's inequality with exponent 
p/{p — 1) and conjugate exponent p that 

2{A + a\\Xt\\i)Z'r^ = 2aZf + 2{A - aC)Zr^ 

p p 

for some e > 0. As for the second term in (4.16) we first estimate 

(r + 7ll^tll?)<7((r/7)'/' + II^J 

Next by the arithmetic-geometric mean inequality we get 

'(r/7)i/2 + ||Xi|ii)'(Zi + ||i^N|U)?'-2 



^ 2 
1 , 



< (^) + ^(in^Niic. + c) + = zr, 

provided that we choose C as the solution to the equation 
'^'^ ' ^ '(||1^N|U + C) = C. 



P V7/ V 
Taken together we thus have 



< I -K„ + I £^2^f + constant, 



where the constant depends on e. Since > we may pick a small enough e such 
that the bracketed expression remains negative. By Gronwall's inequality this then 
proves the result. □ 

For technical reasons, two more a priori moment estimates are required as com- 
plements to Theorems 4.5 and 4.9. 

Lemma 4.10. Let Assumption 4-1 (i) o,nd (iii) hold true. For p > 1, put 

(4.17) MfiXo):^E[\\Xt\h-\\Xo\hr. 
and 

(4.18) MP{Xo):=E\\Xt-Xor. 
Then both Mf and are O (t). 

Proof. From Theorem 4.5 {p > 2) and Propositions 4.2 and 4.3 {p — 1 and 2, 
respectively), we have the estimate 

(4.19) E\\Xt\\{ < WXoir, + W,\\Xo\\{ + B,) 

with explicit expressions for the constants {Pp,Bp\ and with the possibility that 
/3p < whenever the condition in Theorem 4.9 applies. We note a technical point 
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here: this estimate is one-sided in \\Xt\\i — ||^o|li it relies on the signed Assump- 
tion 4-1 (i)- To get a two-sided version we proceed by integrating (2.10), 



t R 



ds 



< 



Cp + CpE\\X,\\P+^ ds ^ O (t) 



by the mean- value theorem, Assumption 4.1 (iii), and various simple estimates 
which we omit for brevity. Using the equivalence of norms, a very similar proof is 
readily obtained for Afl* . □ 

Theorem 4.11 (Initial data perturbation). Let Assumption 4-1 hold true 

and let two solutions of (2.19) Xt and Yt be given, both generated by the same 
counting measure but with different initial data. Put Sq ■= W^o + ^olli- Then 

[Xo ^ Yo] 



{E\\Xt - Yt 



(4.20) 

in terms of 
(4.21) 



||2)i/2 < exp((A/ + AiSo)t) X (^\\Xo - Foil + 
= \\Xo- Foil [1 + {M + /iEo)t + O (t^)] 



-R 



(L' -t- A'So)i + O f<3/2 



R = 



L' + A'So -f O ( 



-(M+t,T.o)s 



and where L' = ||1 N ||oo-^ o.nd similarly for A', and where the indicator function 
[f] = I iff f is true and otherwise. 

Proof. Our starting point is the following differential version of (4.12), 



dt 

(4.22) 

in terms of 



EWXt-Ytf <E 2{M + ^x^o)\\Xt-Yt\ 



[{L' + A'Eo) + l^Ri + A'i?2] ll^t - Yt\ 



Ri = [\\Xt + Yt\U-\\Xo + Yo\\i]\\Xt-Yt\\, 
R2^\\Xt + Yt\\,-\\Xo + Yo\\,. 

Using the properties of the 1-norm over we readily get 

ERl = E (llXtlli - llXolli + ||yt||i - llrolli)' < 2MiiXo) + 2M^{Yo) ^ O (t) , 

^ ^ ^ 

= :-R3 

by Lemma 4.10. The term Ri needs slightly more work. Subtracting and using the 
previous bound we get 



ERl<R3\\Xo-Yo\\^ 
(4.23) > 



E{2[\\Xt\\i-\\Xo\h]^ + 2[\\Yt\\,-\\Yo\ 



1/2 



{E[\\Xt~Yt\ 



iXo-Yofrf' 
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by the Cauchy- Schwartz inequahty. For brevity, put Xt := Xt ~ Xq and Yt := 
Yt — lo- By the properties of the inner product we have that 

\\\Xt - Ytf - ||Xo - Foin ^\\\Xt- Ytf + 2{Xt - Yt,X^ - Yo)\ 

< \\Xt-Yt\\^ + 2\\Xt-Yt\\\\Xo-Yo\\ 

< 2(1 + \\X„- YoimXtW + 2(1 + \\Xo- YoWmr, 
using the integer inequahty. Inserting this into (4.23) and simphfying a bit we get 

ERl < i?3 11^0 - Yof + 8{E[\\Xth - \\Xohf + - rolli]')'^' 

X {l + \\Xo-Yo\\){E\\Xtf + E\\Yt\\Y' 

= i?3 11^0 - Yof + 8(1 + \\Xo- Yo\\)[Mf{Xo) + Mf (Fo)]'/' 

X [M^iXo) + M^{Yo)]'/' < (1 + ||Xo - Foil') O (t) , 

where Lemma 4.10 was used. Starting from (4.22), using the Cauchy-Schwartz 
inequahty several times and combining the bounds for i?i and R2 we infer that 

^E\\Xt - < 2(M + fiJ:o)E\\Xt ~Ytf 
at 

+ [Xo ^ Yo] [l' + A'Eo + + \')0 ^E\\Xt - Yt\\Y^\ 

For the lower order terms we have used the observation that by path-wise unique- 
ness, we may introduce the indicator function according to E\\Xt — Yt\\ = [Xq ^ 
Fol^^ll^t - Yt\\ < [Xq ^ YQ\{E\\Xt - rt||2)i/2. Without this variant of the integer 
inequahty there would be a growing deviation between the two trajectories even 

with identical initial data. The theorem now follows as an application of Lemma 4.8 
with u(t)2 = £;||Xt - Yt||2_ □ 

We next consider perturbations in the reaction coefficients. Given the linear 
dependence on the coefhcients fc^, r — 1...4 in the elementary reactions (2.4)- 
(2.7) a suitable model seems to be that a perturbation kr ^ kr + Sj- in a propensity 
Wr{x) spreads linearly in a relative sense, 

(4.24) \wr{x,kr) — Wr{x,kr + Sr)\ < Constant x SrWr{x,kr). 
We make this formal by simply requiring that 

(4.25) \\w{x)-wsix)h<S\\wix)\U, 

where (5 is a suitable measure of the total perturbation vector and where the per- 
turbed propensity vector function is given by 

(4.26) wsix) = [wi{x, ki + Si), . . . , w_r(x, ku + SR)f. 

We conveniently assume that the entire statement of Assumption 4.1 carries over to 
the perturbed system, and for convenience we shall also assume that all constants 
are the same. By the triangle inequality and Assumption 4.1 (ii) and (iii) we have 
the bound 

\\w{x) - ws{y)\\i < \\w{x) ~ w{y)\\i + \\w{y) ~ ws{y)\\i 

(4.27) < iL + X\\x + y\\i)\\x-y\\+SiT + j\\y\\l). 
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Also, induced from (4.25), it will be convenient to define Sp to be the smallest 
constant such that 



(4.28) 



\\F{x) - Fs{x)\\ = \\-N[w{x)-W5ix)]\\ < Sf\\w{x)\\ 



We shall need a slight extension of Lemma 4.8 as follows. Consider the differential 
inequality d/dtu{t)'^ < a{t)u+ /Su^ +^(t) with a{t) and j{t) non-negative functions, 
but again /3 a signed constant. Using the integrating factor we obtain 

(4.29) 

(e-''/2*u(t))' < u(0)2 + ^ e-^'j{s)ds + e-^/^'a{s) (^e-^/^'u{s)'^ ds, 
written such that the integral version of Lemma 4.8 applies. The final result is 

rt \ 1/2 , „t 



(4.30) u{t) < e^/2t 



We also introduce a convenient notation to keep track of estimates involving 
exponentials and define exp_|_(at) = exp(0 V at) as a majorant to exp(ai). For 
t > and / positive and non-decreasing we get the convenient rule 



(4.31) 



/(s) exp_j,(as) ds < /(t) t exp_|_(Q:i). 



This notation behaves intuitively in most cases, a notable exception being when 
mixing it with ordinary exponentials, e.g., exp(at) exp_,_(— at) — exp_(.(at). Another 
example is that (4.30) takes the compact form 



(4.32) 



u{t) < exp+(/3/2t) 



provided that 7 and a are non-decreasing and that u{0) = 0. This form is used in 
the proof below. 

Theorem 4.12 {Coefficient perturbation). Let two solutions of (2.19) Xt and Yt 
be given, both generated by the same counting measure and with the same initial 
data, but with the propensities for Yt perturbed by S as indicated in (4.25) -(4.26). 
Put Wo := F -|- 7II A"o||x- Then, firstly, there is the "small perturbation", or "short 
time" estimate 

{E\\Xt - YtWy/^ < exp+((Af + L'/2 + {2^i + X')\\Xo\\i)t 



(4.33) 



iS'Wo)'/h'^^ + SFWot + o{t'/^) 



When t is large a generally sharper estimate is 

{E\\Xt~Yt\\Y^^ <e^p+({M + 2fi\\Xo\\i)t 
(4.34) 



where L' and X' are defined in Theorem 4-^1, ci'n-d similarly, 5' — ||1 N ||oo'5- 

The first estimate (4.33) is theoretically the most important one since it shows 
continuity as (5 — ?> 0. 



(5'VFo)i/2^i/2 ^ ^^^^^ ^ + A'||Ao||i)t + O (t"" 
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Proof. We start from a differential version of (4.11), but taking the perturbation 
of the coefficients into account, 

d 



dt " 
(4.35) 



E 



< E 



2{Xt~YuF{Xt)~Fs{Yt)) 
+ {l^n')\w{Xt)-ws{Yt)\ 
2{M + fi\\Xt + Yt\\,)\\Xt ^Ytf + 2SF\\Xt - Yt\\{T + ^\\Yt\\l) 
+ [L' + A'liXt + Yt\\i)\\Xt -Yt\\+ 5'{T + i\\Yt\\\) 



after using the Cauchy- Schwartz inequahty and the bounds (4.27) and (4.28), re- 
spectively. In the proof of Theorem 4.11 we already have constructed bounds for 
two of the terms in the above expression. For the other two terms we use firstly 
from Propositions 4.3 and 4.4 that 

E{T + -f\\Yt\\l)<Wo + 0[t). 

Secondly, writing 

Emwl - \\x,\\\f = Emwt \\x,\\t) 2\\x,\\,E{m\i \\x,\\i) = o{t), 

by Lemma 4.10. Using this and the Cauchy-Schwartz inequality we see that 

E[\\Xt-Yt\\{T + -^\\Y4l)\ < [Wo+o{t'/^)]iE\\Xt-Ytfy/^ 

The final result is 
d 



dt 



E\\Xt - Ftlr < S'Wo + S'O (t) + 2{M + 2^i\\Xo\\l)E\\Xt - FJ' 



L 



' + 2A'||Xo||i + 2SfWo + (/i + A' + 5f)0 (t^/^)] iE\\Xt - Ytf)'/^. 



The asymptotically sharper estimate (4.34) now follows as an application of the 
extended Gronwall's inequality as expressed in (4.32). The short time estimate 
(4.33) is a result of using the integer inequality 

{L' + 2A'||Xo||i) E\\X, - nil < [L' + 2A'||Xo||i) E\\X, - Y,f, 

rather than the Cauchy-Schwartz inequality, and then applying (4.32). This has 
the effect of transferring all terms which do not go to zero with the perturbation S 
into the exponent and is therefore sharper when S and/or t are small. □ 

4.5. Application: deterministic vs. stochastic sensitivity. Inspecting the 
proof of Theorem 4.12, one sees that the essential difference to the corresponding 
deterministic setting is the appearance of the (non-signed) Lipschitz-constants L 
and A. These enter implicitly in (4.35) as the final term and can be traced back 
to Dynkin's formula (2.23). It is therefore interesting to investigate if one can 
construct a simple model displaying a significant sensitivity difference between the 
stochastic and the deterministic formulations ((2.19) and (2.3), respectively). 

This is a nonlinear property, and we see from Table 3.2 that the most expansive 
reaction incorporating only two species is 



(4.36) 



C 



E E, 



where we may think of E as an enzyme and C an intermediate complex. 
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We combine (4.36) with 




(4.37) 



and enforce (C) := E[Coo\ = 10 = -E[-Eoo] {E) at steady-state (note that {E) = 
CiE/ Pe)- We choose units of time such that /Sp = 1 and use the same value of Pe- 
We also conveniently determine ac from the approximation (C) ~ ac/{l3c + k{E)). 
To close the model we make the choice k = 10^ since a relatively large value will 
increase the sensitivity to parameter perturbations (by Table 3.2 we have < fc/4 
in Assumption 4.1 (iv), implying a large exponent in both (4.33) and (4.34)). 

The chemical network thus constructed in a 'bottom-up' fashion is very nearly 
the same as the one considered in [24]. The main difference is the interpretation 
of (4.36), which in [24] is understood to be inhibitive; hence the final product is 
obtained from the intermediate complex by maturation. 

Consider now the effect in C of a perturbation in the amount of enzyme E. 
Specifically, we perturb the system according to aE Qf£;(l — S) otEj"^ such 
that by linearity of (4.37), {E) gets reduced by a factor of 1/2. In an attempt to 
estimate the effect of this perturbation, one might consider to effectively linearize 
the model by replacing (4.36) with 



thereby also making the dynamics one-dimensional. In Figure 4.1 we investigate 
the difference between the nonlinear and linear models by solving numerically the 
deterministic rate equation (2.3) (which is an exact equation in the linear case). 
Although the latter model reacts quicker, the final response is in both cases to 
increase the amount of complex C by a factor of 2, which is what is also intuitively 
expected. 



(4.38) 



C 



kc-{E) 



20 - 




15 - 



O 



10 



5 - 







2 



3 



4 



5 



Figure 4.1. Solid: response for the rate-ODE (2.3) applied to 

(4.36) -(4.37). Dashed: the response for the linear simplified model 

(4.37) -(4.38) is much faster, but both models settle for the same 
asymptotic response. 



However, one can show using the techniques in Propositions 3.5 and 3.6 that the 
original system (4.36)-(4.37) satisfies the bounds {M,y) < (1,25), while the hnear 
version (4.37)-(4.38) admits (M,/i) < (-1001,0) in Assumption 4.1 (iv). Although 
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in the former case the actual figures might be overly pessimistic, the contractive 
character of the latter clearly hints that the linearized model cannot be used to 
predict the sensitivity of the original model. This is confirmed in Figure 4.2 where 
the result of stochastic simulations of (4.36)-(4.37) is displayed. Not only is the 
initial response faster than observed in the nonlinear rate equations in Figure 4.1, 
the final response is also much larger (by a factor of two). The results are even 
more dramatic when, as in Theorem 4.12, one looks at the RMS difference between 
C{t) and the perturbed trajectory Cs{t) (cf. Figure 4.3). 



40 
30 

O 

20 

10( 




t 



Figure 4.2. Circles: response to perturbation of the original non- 
linear model, dotted: unperturbed version. Contrary to the esti- 
mates in Figure 4.1 (which are reproduced here), a fourfold in- 
crease of the amount of complex C is observed. Sample average of 
Af = 10** simulations. 




Figure 4.3. Circles: sample estimate for the RMS difference be- 
tween C{t) and the perturbed trajectory Cs{t). Dotted: ±1/-\/M 
one standard deviation for M = 10"* samples. 
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Although for this example the estimates in Theorem 4.12 yield a qualitative 
insight, they are also very crude compared to the observed results. Since we have 
already commented that the corresponding estimates in the deterministic case are 
quite similar under the present assumptions, this holds true in the plain ODE 
setting as well. To improve on this 'worst case' character of the estimates it seems 
difficult to do more than careful studies on a case-by-case basis. 

5. Conclusions 

We have proposed a theoretical framework consisting of a priori assumptions 
and estimates for problems in stochastic chemical kinetics. The assumptions are 
strong enough to guarantee well-posedness for a large and physically relevant class 
of problems. Explicit upper bounds for the effect of perturbations in initial data 
and in rate constants have been obtained. The assumptions are constructive in 
the sense that explicit techniques for obtaining all postulated constants have either 
been worked out in detail or at least indicated. 

In the course of motivating our setup we have seen that most problems do not 
admit global Lipschitz constants and that one-sided versions do not provide a better 
alternative. Another conclusion worth highlighting is that it pays off to consider 
jump SDEs in a fully discrete setting in that there are potential complications 
in proving path- wise uniqueness when the state space is continuous. A practical 
implication is that some care should be exercised when forming continuous approx- 
imations to these types of jump SDEs. 

It should be clear that we have not considered all possible perturbation results. 
For instance, bounds of the form Esupg^^ || • || were restricted to mere existence in 
the 2nd order moment. Also, ergodicity has not been discussed at all although we 
were able to show under certain conditions that the systems admit asymptotically 
bounded moments. 

For future work we intend to re- visit certain classical results from the perspective 
of the framework developed herein; for example, thermodynamic limit results, time 
discretization strategies, and quasi-steady state approximations — all of which have 
a practical impact in a range of applications. 
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